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Abstract. In some models of quantum gravity, space-time is thought to have a foamy structure with non-trivial optical prop- 
erties. We probe the possibility that photons propagating in vacuum may exhibit a non-trivial refractive index, by analyzing 
the times of flight of radiation from gamma-ray bursters (GRBs) with known redshifts. We use a wavelet shrinkage procedure 
for noise removal and a wavelet 'zoom' technique to define with high accuracy the timings of sharp transitions in GRB light 
curves, thereby optimizing the sensitivity of experimental probes of any energy dependence of the velocity of light. We apply 
these wavelet techniques to 64 ms and TTE data from BATSE, and also to OSSE data. A search for time lags between sharp 
transients in GRB light curves in different energy bands yields the lower limit M > 6.9 • 10 15 GeV on the quantum-gravity scale 
in any model with a linear dependence of the velocity of light oc E/M. We also present a limit on any quadratic dependence. 
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astro-ph/02 10124 complete mathematical model for quantum gravity, and there 
ACT-02-08, CERN-TH/2002-258, MIFP-02-03 are many different approaches to the modelling of space-time 
$— i ' October 2002 foam. Several of these approaches suggest that the vacuum ac- 



1. Introduction 



quires non-trivial optical properties, because of gravitational 
recoil effects induced by the motion of energetic particles. 
In particular, it has been suggested that these may induce a 
In standard relativistic quantum field theory, space-time is con- non-trivial refractive index, with photons of different energies 
sidered as a fixed arena in which physical processes take place, travelling at different velocities. Such an apparent violation 
The characteristics of the propagation of light are considered of Lorentz invariance can be explored by studying the prop- 
as a classical input to the theory. In particular, the special and a g ation of particles through the vacuum, in particular photons 
general theories of relativity postulate a single universal veloc- emitted by distant astrophysical sources (Amelino-Camelia et 
ity of light c. However, starting in the early 1960s (Wheller al - 1998 )- In some quantum-gravity models, light propagation 
1 963), efforts to find a synthesis of general relativity and quan- ma Y also de P end on the P hoton polarization (Gambini & Pullin 
turn mechanics, called quantum gravity, have suggested a need 19 ")> inducing birefringence. Stochastic effects are also pos- 
for greater sophistication in discussing the propagation of light sible > S ivin g rise t0 an energy-dependent diffusive spread in 
in vacuum m e velocities of different photons with the same energy (Ford 

A satisfactory theory of quantum gravity is likely to require 1995, Ellis et al. 2000a). 



a drastic modification of our deterministic representation of 
space-time, endowing it with structure on characteristic scales 
approaching the Planck length £p m p l . There is at present no 



One may discuss the effects of space-time foam on the 
phase velocity, group velocity or wave-front velocity of light. 
In this paper, we discuss only the signature of a modification of 
Send offprint requests to: A.S. Sakharov the group velocity, related to a non-trivial refractive index n(E): 
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v(E) — cj n(E). This may be derived theoretically from a (renor- 
malized) effective Maxwell action r e ff[E,B], where E and B 
are the electric and magnetic field strengths of the propagating 
wave, in the background metric induced by the quantum gravity 
model under consideration. Once the effective Maxwell action 
is known, at least in a suitable approximation, one can analyze 
the photon dispersion using the effective Maxwell equations 
(Ellis et al. 2000b). 

One generally considers the propagation of photons with 
energies E much smaller than the mass scale M characterizing 
the quantum gravity model, which may be of the same order as 
the Planck mass M P , or perhaps smaller in models with large 
extra dimensions. In the approximation E <k M, the distortion 
of the standard photon dispersion relation may be represented 
as an expansion in E/M: 



E 2 = + £i(k/M) + &(k/M? + ...), 



(1) 



which implies the following energy dependence of the group 
velocity 



v*c(l -/(£)). 



(2) 



Here the function f(E) indicates the difference of the vacuum 
refractive index from unity: n{E) — 1 + f(E), which is defined 
by the subleading terms in the series Eq. 

Different approaches to the modelling of quantum gravity 
suggest corrections with different powers of E/M. One of the 
better developed is a string-inspired model of quantum space- 
time, in which the corrections in Eq. start with the third 
power of k/M, as suggested by one particular treatment of D 
branes (Ellis et al. 1997; Ellis et al.1998; Ellis et al. 2000c). 
In this approach, the related violation of Lorentz invariance 
is regarded as spontaneous, and is due to the impacts of light 
but energetic closed-string states on massive D(irichlet) parti- 
cles that describe defects in space-time. In the modern view 
of string theory, D particles must be included in the spectrum, 
and hence also their quantum fluctuations should be included 
in a consistent formulation of the ground-state vacuum. In the 
model of Ellis et al. 1997; Ellis et al. 1998; Ellis et al. 2000c, 
the scattering of the closed-string state on the D particle in- 
duces recoil of the the latter, which distorts the surrounding 
space-time in a stochastic manner, reflecting the foamy struc- 
ture of space-time. 

In such a picture, the recoil of the massive space-time def- 
fect, during the scattering with a relativistic low-energy probe 
such as a photon or neutrino, distorts the surrounding space- 
time, inducing an effective net gravitational field of the form 



Gl>i ~ M 



(3) 



The dispersion-relation analysis (Ellis et al. 2000b) of the 
Maxwell equations in the non-trivial background metric per- 
turbed by such a gravitational field results in a linear depen- 
dence of the vacuum refractive index on the energy: 



f(E) 



(-)• 

\M) 



In some other realisations of quantum gravity, odd powers of 
k/M in Eq. Q may be forbidden (e.g. Burgess et al. 2002) 
by selection rules such as rotational invariance in a preferred 
frame. In this case, the leading correction to the refractive index 
takes the form 



f(E) 



\Ml 



(5) 



(4) 



We assume that the prefactors in both cases are positive, reflect- 
ing the fact that there should be no superluminal propagation 
(Moore & Nelson 2001). This requirement is not necessarily 
respected in some models based on the loop-gravity approach 
(Gambini & Pullin 1999a; Alfaro et al. 2000). 

The study of short-duration photon bursts propagating over 
cosmological distances is a most promising way to probe this 
approach to quantum gravity (Amelino-Camelia et al. 1998): 
for a recent review, see (Sarcar 2002). The modification of 
the group velocity Eq. (0 would affect the simultaneity of the 
arrival times of photons with different energies. Thus, given a 
distant, transient source of photons, one could measure the dif- 
ferences in the arrival times of sharp transitions in the signals 
in different energy bands. Several different types of transient 
astrophysical objects can be considered as sources for the pho- 
tons used to probe quantum-gravity corrections such as Eq. (0} 
and Eq. to the vacuum refractive index (Amelino-Camelia 
et al. 1998; Ellis et al. 2000b; Biller et al. 1999; Schafer 1999). 
These include Gamma-Ray Bursters (GRBs), Active Galactic 
Nuclei (AGNs) and pulsars. 

A key issue in such probes is to distinguish the effects of 
the quantum-gravity medium from any intrinsic delay in the 
emission of particles of different energies by the source. Any 
quantum-gravity effect should increase with the redshift of the 
source, whereas source effects would be independent of the 
redshift in the absence of any cosmological evolution effects 
(Ellis et al. 2000b). Therefore, in order to disentangle source 
and propagation effects, it is preferable to use transient sources 
with a known spread in redshifts z. At the moment, one of the 
most model-independent ways to probe the time lags that might 
arise from quantum gravity is to use GRBs with known red- 
shifts, which range up to z ~ 5. 

Increasing numbers of redshifts have been measured in re- 
cent years, and the spectral time lags of GRB light curves 
have been investigated in a number of papers (Norris et al. 
1994; Band 1997; Norris et al. 2000; Ellis et al. 2000b; Norris 
2002). It is important to detect quantitatively temporal struc- 
tures which are identical in different spectral bands, to compare 
their time positions. Unfortunately, pulse fitting is problematic 
(Norris et al. 1996; Ellis et al. 2000b; Norris et al. 2000) in the 
cases of many bursts, because of irregular, overlapping struc- 
tures in the light curves. As a result these studies often lack 
the accuracy to characterize short-time features in the bursts 
that are evident to the eye. The cross-correlation method (Band 
1997; Norris et al. 2000; Norris 2002) does not use a rigorous 
definition of a spike in a pulse; it relies, instead, on a calcula- 
tion of the cross-correlation functions (CCFs) between differ- 
ent spectral bands directly in the time domain. However there 
are some ambiguities in the interpretation of CCF peaks, which 
can lead in some cases to unclear conclusions about the spectral 
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evolution. In particular this is the case when a GRB light curve 
contains an emission cluster of closely spaced spikes (e.g. spac- 
ing of order of the width of a spike); then the width of the CCF's 
central peak, the position of which actually measures the spec- 
tral time lag (Band 1987), may reflect the duration of the whole 
cluster and not of the individual spikes, whereas only the nar- 
row individual constituents (spikes) of such an emission cluster 
can mark with a good accuracy the arrival time of radiation, so 
as to apply for a search for quantum gravity effects (Amelino- 
Cameliaetal. 1998). 

In this paper, we seek to overcome the problems mentioned 
above by using wavelet transforms to remove noise, to resolve 
overlapping structures and to classify quantitatively the irregu- 
larities of GRB light curves with known redshifts. The ability 
of the wavelet technique to characterize burst morphology al- 
lows us to improve significantly the accuracy of the measure- 
ments of time lags, independently of the degree of spikes sepa- 
ration inside the emission clusters, 

increasing the sensitivity to quantum-gravitational 
corrections. We analyse the light curves of GRBs with 
known redshifts triggered by the Burst And Transient 
Source Experiment (BATSE) aboard the Gamma 
Ray Observatory (GRO) (see the GRO webpage at 
http : //cossc . gsf c . nasa . gov/cgro/index . htmli, 
searching for a redshift dependence of spectral time lags 
between identical sharp signal transitions detected by wavelet 
transforms in different spectral bands. For several GRBs 
among the triggers under consideration, one can compare the 
BATSE light curves with those measured at higher energies 
by the Orientated Scintillation Spectrometer Experiment 
(OSSE) aboard the GRO. We also demonstrate that the wavelet 
technique can deal with the leading parts of the GRB light 
curves recorded by the BATSE time trigger event (TTE), which 
improves the time resolution substantially. Unfortunately, in 
all the cases except GRB980329, the TTE data do not cover 
enough of the light curve to exhibit coherent structures in 
different spectral bands, which would increase the sensitivity 
to higher quantum-gravity scales. 

We find that the combination of all the available data, when 
analyzed using wavelet transforms, prefers marginally a lin- 
ear violation of Lorentz invariance between 10 15 GeV and 
10 16 GeV, although the effect is not significant. We prefer to in- 
terpret the data as giving a limit on the linear quantum-gravity 
scale: 



M > 6.9 ■ 10 n GeV, 



(6) 



which we consider to be the most robust and model- 
independent currently available. 

The layout of this paper is as follows. In Section 2 we dis- 
cuss the propagation of light in an expanding Universe, estab- 
lishing the basic formulae we use subsequently in our analy- 
sis of time lags. The fundamental definitions and features of 
wavelet transforms are reviewed in Section 3, and we describe 
in Section 4 how wavelet shrinkage can be used to remove 
noise from GRB spectra. The 'zooming' technique for localiz- 
ing variation points in GRB light curves is described in Section 
5, and Section 6 uses this technique to analyze time lags. Our 



limits on linear and quadratic quantum-gravity models are ob- 
tained in Section 7, and we discuss our results in Section 8. In 
addition, Appendix A discusses signal threshold estimation in 
the wavelet approach, and Appendix B recalls some aspects of 
the Lipschitz characterization of singularities. 

2. Light Propagation in the Expanding Universe 

Light propagation from remote objects is affected by the ex- 
pansion of the Universe and depends upon the cosmological 
model (Ellis et al. 2000b). Present cosmological data motivate 
the choice of a spatially-flat Universe: Q to tai = + = 1 
with cosmological constant £2a =s 0.7: see (Bahcall et al. 1999) 
and references therein. The corresponding differential relation 
between time and redshift is 
dz 

(7) 



(l+z)h(z) 
where 



h{z) = Vn A + n M (l+z) 3 . 



(8) 



Thus, a particle with velocity u travels an elementary distance 

udz 



udt = —H n 



(9) 



(l+z)h(z)' 

giving the following difference in distances covered by two par 
tides with velocities differing by Am: 



AL 



- H 1 f A 
- H ° J (1 + 



Audz 



z)h(z) 



(10) 



We consider two photons traveling with velocities very close to 
c, whose present day energies are E\ and Ei. At earlier epochs, 
their energies would have been blueshifted by a factor 1 + z. 
Defining AE = Ei - E\, we infer from equation Eq. (0 that 



Am = 



AE(l +z) 
M 



(11) 



in the case Eq. (|4} of a linear L-dependence of the velocity of 
light, and 

A£ 2 (l +z) 2 



An 



M 2 



(12) 



for the quadratic correction Eq. (|5j- Inserting the last two ex- 
pressions into Eq. dlOi . one finally finds that the induced dif- 
ferences in the arrival times of the two photons with energy 
difference AE are 



At = H, 



and 



,_i A£ C dz 
M J h(zY 



(13) 



(14) 



for the linear and quadratic types of correction, respectively. 

In the following, we look for such time differences in the 
arrival times of photons with energy difference AE propagating 
in such a flat expanding Universe with a cosmological constant 
(BahcaUetal. 1999). 



\ M) J h(z) 
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3. What can Wavelet Transforms do? 

Wavelet transforms (WT) (for a review, see Dremin et al. 
2001) are used to represent signals which require for their 
specification not only a set of typical frequencies (scales), but 
also knowledge of the coordinate neighbourhoods where these 
properties are important. The most important principles distin- 
guishing a wavelet basis from a windowed Fourier transform 
basis are dilatations and translations. Dilatations enable one 
to distinguish the local characteristics of the signal at various 
scales, and translations enable one to cover the whole region 
over which the signal is studied. 

The wavelet transform of a function / at the scale s and po- 
sition u is computed by convoluting / with a wavelet analyzing 
function: 

b 

Wf(u,s) = J" f(t)y us (t)dt. (15) 

a 

The analyzing function ifr us is obtained through dilatation by a 
scale factor s and translation by an amount u from a basic (or 
mother) wavelet if/: 



<M0 = 




It is obvious that if/ must satisfy an admissibility condition 
which guarantees the invertibility of the wavelet transform. In 
most cases, this condition may be reduced to the requirement 
that if/ is a function with zero mean (Mallat 1998): 

oo 

j if/(t)dt = 0. (17) 

-co 

In addition, if/ is often required to have a certain number of 
vanishing moments: 

CO 

j t"if/(t)dt = 0, n = 0,1,...,/?. (18) 

— oo 

In general, this property improves the efficiency of if/ for de- 
tecting features (singularities) in the signal, since it is blind to 
polynomials up to order N. One may say that the action of s on 
the function if/, which must be oscillating according to Eq. d!7i . 
is a dilatation if s > or a contraction if s < 0. In either case, 
the shape of the function is unchanged, it is simply spread out 
or squeezed. 

A transform Eq. (I15> over a suitable wavelet basis is usu- 
ally called a continuous wavelet transform (CWT). A wavelet 
transform whose wavelets if/ are constructed in such a way that 
the dilated and translated family 

1 lt-Vm\ 

« )= V^— ) ; (19) 

where j, m are integers, is an orthonormal basis for all functions 
/ satisfying the condition 

j \f\ 2 (t)dt < +oo. (20) 



This is called a discrete wavelet transform (DWT): for a review, 
see (Dremin et al. 2001; Mallat 1998). 

The CWT is mostly used for the analysis and detection of 
signals, whereas the DWT is more appropriate for data com- 
pression and signal reconstruction. Combining these two types 
of wavelet transforms provides an advanced technique for pick- 
ing up the positions of particular breaks in the structures of ob- 
served GRB light curves in different energy bands, which we 
use here to look for non-trivial medium effects on the propaga- 
tion of photons due to quantum gravity. 

Orthogonal wavelets Eq. dl9> dilated by factors V are sen- 
sitive to signal variations with resolutions 7rK This property 
can be used to make a sequence of approximations to a signal 
with improving resolutions (e.g. Mallat 1998). For a function 
satisfying the condition Eq. j20L the partial sum of wavelet 
coefficients YX=-°a dj,k^j,k can be interpreted as the difference 
between two approximations to / with resolutions 2~- ,+1 and 
2~ J . Adapting the signal resolution allows one to process only 
the details relevant to a particular task, namely to estimate in- 
tensity profiles of GRB light curves preserving the positions of 
sharp signal transients. 

CWTs can detect with very high precision the positions 
where the intensity profile of a GRB light curve, as estimated 
by a DWT, changes its degree of regularity. Since if/ has zero av- 
erage, a wavelet coefficient Wf(u, s) measures the variation of 
/ in a neighborhood of u whose size is proportional to s. Sharp 
signal transitions create large-amplitude wavelet coefficients. 
As we see in the following section, the pointwise regularity of 
/ is related to the asymptotic decay of the wavelet transform 
Wf(u, s) when s goes to zero. Singularities are detected by fol- 
lowing across different scales the local maxima of the wavelet 
transform. We use this 'zooming' capability to define the po- 
sitions of mathematically similar transients (irregularities) in 
GRB light curves observed in different energy bands. These 
therefore provide the best information about the arrival times 
of photons associated with universal intrinsic emission features 
at the sources. 

4. Extraction of the GRB Intensity Profiles by 
Wavelet Shrinkage 

The observed GRB light curves typically feature a relatively 
homogeneous, nonzero background intensity, above which 
some inhomogeneous structure is apparent (Kolaczyk 1997). 
In the following, we demonstrate that when such a temporally 
inhomogeneous signal as the light curve of a GRB contains 
both structure and noise, the ability of the DWT to compress the 
information in this signal leads efficiently to a simple but effec- 
tive noise removal procedure. This wavelet shrinkage technique 
(Donoho 1993; Donoho et al. 1995), based on the thresholding 
of the DWT, allows one to separate the structure of the signal 
from the noise, whilst retaining information about the position 
of irregularities of the signal, as provided by the support of the 
mother wavelets. 

In practice, DWTs break a function down into a coarse ap- 
proximation at a given scale, that can be extended to succes- 
sive levels of residual detail on finer and finer scales. The full 
decomposition may be expressed in terms of the scale func- 
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Fig. 1. Data from GRB990308 with z = 1 -60 collected by BATSE trigger 7457 in a total number of 2 11 64 ms bins in the energy 
band 115-320 keV (top panel), and its Symmlet-10 (Mallat 1998) discrete wavelet transform (DWT) at the level L = 6. The 
horizontal axis corresponds to the same time scale as in the original burst. Each level on the vertical axis shows the wavelet 
coefficients at a given resolution level (i.e., scale). The wavelet coefficients are represented by spikes whose size and direction 
(up or down) is determined by the magnitude and sign (+ or -) of the coefficient. 



tion <f>j JX and the discrete wavelet i//p, discussed already. The 
scale function looks much like a kernel function, and a finite 
linear combination of dyadic shifts of this function provides a 
generic coarse approximation. Further linear combinations of 
dyadic shifts of the wavelet function supply the residual detail. 
By considering sequences of DWTs with increasing numbers 
of dyadic dilatations, the detail at each of the corresponding 
successive scales is recovered. 

We represent the light curve of a given GRB by a binned 
discrete signal {Xo,Xu . . . ,X„-{\ of dyadic length n = 2 J , 
where J > 0. The DWT of such a signal results in a vector of 
length n of wavelet coefficients. The signal is said to have been 
sampled at level J. At some coarser resolution level (scale) 
L < J, the wavelet coefficient vector contains 2 L scale coef- 
ficients Clj), . . . , c L 2 l-i and 2 J detail coefficients djjo, . . . , dj^j-\ 
at each of the levels j = L,...,J— 1. Fig.[2displays a typical 
example with 7=11 and L — 6. 

Observations of a given GRB light curve can be represented 
by the sum 



where the intensity profile f[n] is contaminated by the addition 
of noise, which is modelled as a realization W[n] of a random 
process whose probability distribution is known. The intensity 
profile / is estimated by transforming the noisy data X[n] with 
the 'decision operator' D: 



F = DX. 



(22) 



X[n] = f[n] + W[h\, 



(21) 



The 'risk' of the estimator F off is the average loss, calculated 
with respect to the probability distribution of the noise. The 
numerical value of the risk is often specified by the signal-to- 
noise ratio (SNR), which is measured in decibels. 

To reduce the noise level of W, while preserving the degree 
of regularity of the intensity profile /, we use a soft thresh- 
olding procedure. This procedure sets to zero all coefficients 
smaller in magnitude than some threshold T, and shrinks coeffi- 
cients larger than T towards zero by amounts T, as described in 
more detail in AppendixlAl This performs an adaptive smooth- 
ing that depends on the regularity of the signal /. In a wavelet 
basis 1 where large-amplitude coefficients correspond to tran- 

1 In general, the thresholding procedure can be applied to any basis 
for the signal decomposition (e.g. Mallat 1998). 
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Fig. 2. The GRB990308 light curve after preprocessing which sets the median absolute deviation of wavelet coefficients at the 
fine scale equal to unity, as described in more detail in Appendix |A] The intensity profile estimated by the wavelet shrinkage 
procedure at the level L = 6 is presented in the bottom panel. 



sient signal variations, this means that the estimator discussed 
in Appendix 1X1 keens only transients coming from the original 
signal, without adding others due to the noise. After the pre- 
processing, which sets the median value of wavelet coefficients 
of the signal at the finest scale to unity, the threshold is esti- 
mated to be T — y2 log n. The example of an intensity profile 
estimated by this wavelet shrinkage procedure, as described in 
AppendixIXl is shown in Fig. |3] 



In general, the wavelet shrinkage procedure Eq. JA.2> de- 
scribed above guarantees with high probability that \d^ m \ < 

\dj m \ (e.g. Mallat 1998), implying that the estimator F is at 
least as regular as the 'original' intensity profile /, because its 
wavelet coefficients have smaller amplitudes. Thus we use this 
property of the DWT of separating very effectively the struc- 
tures in the GRB intensity profiles from noise, in the form 
of two subsets of wavelet coefficients, large and small ones. 
The thresholding procedure deletes wavelet coefficients below 
the threshold value, and diminishes the others by the threshold 
value. This tends to preserve both broad and narrow features, 
while significantly reducing noise fluctuations, after the recon- 
struction of the intensity profile by the inverse DWT. 



5. Detection of Variation Points in GRB Light 
Curves 

A wavelet transform can focus on localized signal structure via 
a 'zooming' procedure that reduces progressively the scale pa- 
rameter. To identify the variability points of two different light 
curves and characterize their structures, it is necessary to quan- 
tify precisely the local regularity of the function which repre- 
sents the intensity profile of the original signals. The appro- 
priate tools are the Lipschitz exponents 2 . These are defined in 
Appendix|B] and can provide uniform regularity measurements 
of a function / not only over time intervals, but also at any point 
v. 

We use the Lipschitz exponent a to characterize variation 
points of the reconstructed intensity profiles of GRBs accumu- 
lated in different energy bands. The comparison of the positions 
of variation points with the same values of a gives the infor- 
mation about the arrival times of photon probes with different 
energies, enabling one to probe for quantum-gravity time-delay 
phenomena. 



2 Lipschitz exponents are also called Holder exponents in the math- 
ematical literature (Dremin et al. 2001; Mallat 1998). 
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Fig. 3. The CWT image (middle panel) of the GRB990308 intensity profile (top panel). The horizontal and vertical axes give, re- 
spectively, the position u and log 2 s (where s is the time scale in seconds), The shadings (colours from yellow to blue) correspond 
to negative, zero and positive wavelet transforms respectively. Singularities create large amplitude coefficients in their cone of 
influence. The modulus maxima (bottom panel) of Wf(u, s) obtained from the matrix of CWT (middle panel) pointing towards 
the time positions of singularities at the fine scale. 



The decay of the CWT amplitude as a function of scale 
is related to the uniform and pointwise Lipschitz regularity of 
the signal. Thus, measuring this asymptotic decay is equiva- 
lent to 'zooming' into signal structures with a scale that goes to 
zero. Namely, when the scale s decreases, the CWTs Wf(u, s) 
Eq. G3 measures fine-scale variations in the neighborhood of 
u. One can prove (e.g. Mallat 1998) that \Wf(u, s)\ decays 
like s a+l l 2 over intervals where / is uniformly Lipschitz a. 
Furthermore, the decay of \Wf(u, s)\ can be controlled from its 
local maxima values. 

A 'modulus maximum' (e.g. Mallat 1998) is any point 
(uo, so) such that \ Wf(u, sq)\ is locally maximal at « = uq. This 
local maximum should be a strict local maximum in either the 
right or the left neighborhood of uq. Any connected curve s(u) 
in the scale-time plane (u, s) along which all points are modu- 
lus maxima, as illustrated in Fig. [3] is called a 'maxima line'. 
Singularities of a function / are detected by finding the abscissa 
where the wavelet modulus maxima converge on fine scales 
(e.g. Mallat 1998). Only at such points can / be singular, i.e., 
with exponent a < 1 . This result guarantees that all singular- 



ities are detected by following the wavelet transform modulus 
maxima at fine scales. Fig. shows an example where all the 
significant singularities are located by following the maxima 
lines. The positions of these singularities are located by the 
modulus maxima lines at the fine scale of decomposition. 

To be sensitive to both sharp and smooth singularities, one 
has to use wavelets with two vanishing moments, so as to gen- 
erate the CWT Eq. (1151 of the reconstructed intensity profiles. 
The most suitable one is the second derivative of a Gaussian 
(Mexican hat) mother wavelet: 



¥(*) = 



f^- llcxpl-— T |. 



2cr 2 



(23) 



because of the property that the modulus maxima of Wf(u, s) 
with the wavelet Eq. d23i belong to connected curves that are 
not broken as the scale s decreases (e.g. Mallat 1998), which 
guarantees that all maxima lines propagate to the finest scales. 
The dilatation step s is generally set to s = 2^ A , where A is the 
number of intermediate scales (voices) for each octave. Thus, 
if the voice lattice is sufficientely fine, one can build maxima 
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Fig. 4. The light curve for GRB990308 obtained by BATSE with trigger 7457, binned with 64 ms resolution in four spectral 
bands. 



lines with very high precision. Connecting maxima into lines as 
in Fig.[3]is a procedure for removing spurious modulus maxima 
created by numerical errors in regions where the CWT is close 
to zero. 

Sometimes the CWT may have a sequence of local max- 
ima that converge to a point v on the abscissa, even if / is per- 
fectly regular at v. Thus, to detect singularities it is not suffi- 
cient merely to follow the wavelet modulus maxima across the 
scales. One must also calculate the Lipschitz regularity from 
the decay of the modulus maxima amplitude. If, for some scale 
s < so, all the modulus maxima that converge to v are included 
in a cone: 



\u — v\ < Cs, 



(24) 



then / has an isolated singularity at v. Conversely, the absence 
of maxima below the cone of influence Eq. ( 1241 implies that / 
is uniform in the neighborhood of any point t + v beyond the 
scale sq. 

The Lipschitz regularity at v is given by the slope 
log 2 \Wf(u, s)\ as a function of log 2 s along the maxima lines 
converging to v, namely 



log 2 |W/(«,s)| 



a + — I log-, s + const. 



(25) 



Actually, the Lipschitz property Eq. (IB. 1 1 approximates a func- 
tion with a polynomial p v in the neighborhood of the point v. 
The CWT estimates the Lipschitz exponents of the function by 
ignoring the polynomial p v itself. Moreover, if the scale sq is 
smaller then the distance between two consecutive singulari- 
ties, to avoid having other singularities influence the value of 
Wf{u, s), and the estimated Lipschitz exponent a + 1/2 < 1.5, 
the function / exhibits a break at v, which can be detected by 
following the modulus maxima chain. 

In this paper, to define significant points in the time se- 
ries of the signal, we do not apply fit functions that select only 
prominent peaks, as was done in (Ellis et al. 2000b; Norris et 
al. 1995). We consider a more general class of relatively sharp 
signal transitions, marked by Lipschitz irregularities picked out 
by the CWT 'zoom' technique, which we call genuine varia- 
tion points. 

6. Analysis of Time Lags in Emissions from GRBs 
with Known Redshifts 

Once the BeppoSAX satellite began to localize long bursts 
in the sky to within a few arcminutes, and distribute their 
locations to observers within hours, it turned to be possible 
to discover X-ray, optical, and radio afterglows (Costa et al. 
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Fig. 5. The estimated intensity profiles of GRB990308 - see Fig. [4]- obtained in four spectral bands by using a Symmlet-10 
(Mallat 1998) basis at level L = 6. The signal-to-noise ratios are at the levels SNRl = 23.07, SNR2 = 22.71, SNR3 = 23.58 and 
SNR4 = 23.87 for each band, respectively. All variation points founded by CWT zoom are marked by circles. The behaviours 
of the Lipschitz exponents a are estimated. Seven pairs of genuine variation points in the first and third spectral bands have been 
detected. 



1997; van Paradijs et al. 1997; Frail et al. 1997), and host 
galaxies. Subsequent observations led to the spectroscopic de- 
termination of GRB redshifts, using absorption lines in the 
spectra of the afterglows and emission lines in the spec- 
tra of the host galaxies. By now, redshifts have been mea- 
sured for about 20 bursts (see for example (Norris et al. 
2000; http://www.aip.de/~jcg/grbrsh. html; Amati et 
al. 2002) and references therein). 

Our first aim is to measure the timings of genuine varia- 
tion points, characterized by Lipschitz exponents as discussed 
above, for different spectral bands in the light curves of dis- 
tant GRBs. Correlating the time lags between different energies 
with the GRB redshifts, we then try to extract time delays re- 
lated to the refractive index that may be induced by quantum 
gravity. 

We use genuine variation points with the same Lipschitz 
exponents a, measured in different energy bands, and assume 
that any initial relative time lags attributable to the properties of 
source are independent of redshift. Thus, the key to disentan- 
gling quantum-gravity effects is reduced to the problem of de- 
tecting genuine variation points with the highest possible pre- 



cision. The biggest uncertainties in our analysis come from our 
procedure for estimating the DWT intensity profiles, whilst the 
errors generated by the CWT zoom are negligible. The error 
in the wavelet-shrinkage procedure is defined by the time-bin 
resolution in the analysis of the light curve and the support of 
the DWT mother wavelet. 

As shown in Table ^ we use GRB light curves from 
BATSE, which have been recorded with 64 ms temporal reso- 
lution in four spectral channels. Unfortunately, the BATSE cat- 
alog of light curves includes only about a half of the GRBs 
with known redshifts. The light curves of the other GRBs 
with known redshifts have been collected by other satellites 
(BeppoSAX, HETE), and the data are not available publicly. 
The BATSE lower-level discriminator edges define the chan- 
nel boundaries at approximately 25, 55, 115 and 320 keV. We 
look for the spectral time lags of the light curves recorded in 
the 1 15 - 320 keV energy band relative to those in the lowest 
25 - 55 keV energy band, providing a maximal lever arm be- 
tween photon energies. We do not use the fourth BATSE chan- 
nel with energies between 320 keV and 1 .9 Me V, as these have 
ill-defined energies and poorer statistics - see Fig. |4] Instead, 



10 



J. Ellis et al.: Quantum-Gravity Analysis of Gamma-Ray Bursts using Wavelets 



2.5 



1.5 



0.5 



< o 



-1 



-1.5 



-2.5 



970508 



991216 



990123 



990510 



990308 



980329 



971214 



980703 



980425 



0.2 



0.4 



0.6 



0.8 



1.2 



1.4 



1.6 



1.8 



K, 



Fig. 6. Spectral time lags between the arrival times of pairs of genuine variation points detected in the third and first BATSE 
spectral bands. The analysis has been done for 9 GRB light curves collected with time resolution 64 ms. The solid line shows the 
best linear fit versus K\. 



we compare the rather more energetic light curves accumulated 
by OSSE 3 with the 1 15 - 320 keV BATSE light curves, which 
increases the lever arm for probing photon propagation into the 
MeV range. 

Since we apply the CWT zoom technique to detect the gen- 
uine variation points of the reconstructed intensity profiles, we 
impose some conditions on the choice of shrinking wavelet 4 . 
In general, when choosing the appropriate wavelet basis, one 
has to strike a balance between the degree of regularity of the 
wavelet, the number of its vanishing moments p and the size 
of its support. It is clear that the size of the support defines 
uncertainties in the positions of genuine variation points af- 
ter the reconstruction of intensity profile. This consideration 
motivates using the DWT basis with the most compact sup- 
port for the wavelet shrinkage procedure. On the other hand, 
to preserve maximally the regularity of the original signal, one 



3 We are grateful to M. Strikman for kindly providing us with OSSE 
data. 

4 In most cases, discrete wavelets cannot be represented by an ana- 
lytical expression or by the solution of some differential equation, and 
instead are given numerically as solutions of functional equations (e.g. 
Mallat 1998). 



should use wavelets with a high degree of regularity. In addi- 
tion, one should avoid disturbing significantly the alignment of 
peaks of the original light curves, which motivates using sym- 
metrical discrete wavelets. 

The discrete wavelet that best reconciles the above require- 
ments is that called Symmlet-p (e.g. Mallat 1998). It is the 
most symmetric, regular discrete wavelet with minimum sup- 
port. The number p of vanishing moments defines the size of 
the support, and consequently the errors of the position estima- 
tions \{2p - 1) x bin - size. Moreover the same number p of 
vanishing moments defines the regularity of Symmlet-p. For a 
large number of vanishing moments, the Lipschitz regularity of 
Symmlet-p is 0.275/j (Daubechies 1991). Thus, to have more 
then 2 continuous derivatives, p should exceed 8. 

For the selected GRBs in Table ^ we nave performed 
the wavelet shrinkage procedure using a Symmlet-10 basis 5 . 
At sufficiently high signal-to-noise ratio levels, this procedure 
tends to preserve the regularity of the light curves. In some 
cases, namely for GRB980329 and GRB970508, we applied 



5 The coefficients of the Symmlet filters are tabulated in WAVELAB 
toolbox 

(http://www-stat.stanford.edu/ "wavelab), for example. 
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Fig. 7. The combination of 64 ms BATSE time-lag measurements shown in Fig. [6] with the measurement obtained from the 
BATSE-OSSE comparison and the TTE portion of the GRB980329 light curve, with resolution 2.7 ms. 



the translation-invariant (Mallat 1998) version of the shrink- 
age procedure with reduced threshold. This procedure implies 
averaging estimates produced from the original signal itself 
and from all shifted versions of the signal, and allows one to 
avoid artefacts while preserving the real transient structure. 
Subsequently, we apply the CWT zoom technique for recon- 
structing intensity profiles to identify the arrival times of gen- 
uine variation points and estimate their Lipschitz exponents in 
every spectral band. We consider that a genuine variation point 
has been detected if it has Lipschitz exponent a < 1 . Genuine 
variation points found in the vicinity of each other, but belong- 
ing to two different spectral bands, are considered to have been 
generated at the source if the values of their Lipschitz expo- 
nents are equal to each other. The other variation points with 
a substantially exceeding 1 exhibit only smooth transitions of 
the signal, and do not mark sharp transient time structures. 
We recall that only sharp transient structures are important in 
the search for spectral time lags. For seven GRBs out of nine, 
we detected more than one pair of identical genuine variation 
point per light curve, as seen in Tabled The systematic errors 
(Kolaczyk 1997) were estimated by using different resolution 
levels (L = 6, 5, 4) in the wavelet shrinkage procedure. 



In order to probe the energy dependence of the velocity of 
light that might be induced by quantum gravity, we have com- 
piled the whole available data in Table ^ as functions of the 
variables K\ and K q , defined by the integrals in Eq. Jl 3i and Eq. 
dl4> . respectively. In the case of linear quantum-gravity correc- 
tions, the variable takes the form 



I 



dz 



(26) 



whilst for the quadratic case we use 



I 



(l+z)dz 
h(z) 



(27) 



Since both Eq. dl 31 and Eq. dl4t exhibit linear dependences on 
the variables Eq. J26t i and Eq. J27I respectively, we perform 
a regression analysis for a linear dependence of the time lags 
between pairs of genuine variation points, in the form 



At = aK + b. 



(28) 
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Fig. 8. The x 1 function defined in Eq. J35I as a function of the quantum gravity scale, for the scenario with a linear energy 
dependence of the vacuum refractive index, as obtained using different combinations of data sets. The solid line, which is used to 
establish the lower limit Eq. d36l >. corresponds to the combination of 64 ms BATSE data with both TTE and OSSE-BATSE data. 



The result of our regression fit to the full 64 ms statistics for 
linear quantum-gravity corrections Eq. (0} is shown in Fig. [6] 
The best fit corresponding to Fig.[6]is given by 

At = 0.60(±0.46)tfi - 0.72(+0.53). (29) 

Following the same procedure for the quadratic corrections, 
one gets 

At = 0.17(±0.17)# q - 0.42(±0.32). (30) 

More precise results can be obtained by combining BATSE 
and OSSE data. Four light curves accumulated by OSSE ex- 
ibit structures that can be compared with similar features ob- 
served by BATSE, as seen in Tabled Since the OSSE data are 
at higher energies: 0.15 - 10 MeV, one has to rescale the re- 
sults of OSSE-BATSE comparison in order to combine them 
directly with results obtained using the third BATSE channel, 
by a factor: 
110keV-55keV 

, (31) 

3 MeV- llOkeV 

where 3 MeV is the energy at which the contribution of flux 
accumulated by OSSE becomes significant. The spectral infor- 
mation for GRB980123 indicate that half of the total flux has 
been accumulated in the energy range 3-6 MeV. 



One may also increase the sensitivity of the determination 
of time lags by using BATSE time-tagged event (TTE) data, 
which record the arrival time of each photon with a precision 
of 2 yus, in the same four energy channels. The onboard mem- 
ory was able to record up to 32,768 photons around the time of 
the BATSE trigger. Typically, this quota of photons was filled 
in 1 or 2 s. For short GRBs, the mean structure of the whole 
light curve might be in the TTE data, along with substantial 
periods of background emission after the burst, whilst for the 
long-duration GRBs that we analyze, as in Table ^ the TTE 
data cover only the leading portion of light curve. We have re- 
binned with resolution ^ 1 ms the leading portions of all the 
GRBs from Table^using TTE data. Only one light curve, that 
of GRB980329, yields a signal with clearly identified isolated 
singularities in the first and third spectral bands. The statistics 
available to detect genuine variation points in this light curve 
yield a resolution of 2.7 ms. Combining this TTE measurement 
with the 64 ms BATSE measurements and the results of our 
BATSE-OSSE comparisons, we get the following results: 

Af = 0.010(±0.022)£i - 0.053(+0.026) (32) 
At = 0.003(±0.006)tf q - 0.048(+0.016) (33) 
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Fig. 9. The same^ 2 function as in Fig. [8] which is defined in Eq. (I35L now plotted as a function of the inverse quantum-gravity 
scale. The solid line corresponds the final lower limit Eq. d36l >. 



for linear Fig. [7] and quadratic corrections respectively. The 
single BATSE point with much higher precision than the oth- 
ers does not improve substantially the significance of the fits. 
However, including the OSSE-BATSE and TTE measurements 
into the overall fit does improve the sensitivity (see the next 
section) and makes the result less dependent on the properties 
of individual sources. 

The leading parts of other light curves from Tabled which 
do not exhibit coherently variable structures, can be character- 
ized as fractal signals without isolated singularities. One can 
also analyze such singularities with CWT (e.g. Dremin et al. 
2001; Mallat 1998), but such a study lies beyond the scope of 
this paper. 



7. Compilation of Limits on Quantum Gravity 

We now analyze the likelihood function to derive the results of 
our search for a vacuum refractive index induced by quantum 



gravity. We establish a 95 % confidence-level lower limit on 
the scale M of quantum gravity by solving the equation 

oo 

^ = 0.95 (34) 



where oo symbolizes a reference point fixing the normalization. 
In our case, we choose as reference point M = 10 19 GeV, the 
Planck mass, which corresponds to the highest possible scale 
above which quantum-gravity effects vanish and corrections to 
the vacuum refractive index become infinitesimally small. In 
practice variations of reference point, even by an order of mag- 
nitude, does not influence the final result. 

We use the fact that only the coefficient a in Eq. ( 12 8 1 is 
related to the quantum-gravity scale, whereas b includes a pos- 
sible unknown spectral time lag inherited from the sources, 
which we assume to be universal for our data set. With this 
assumption, one can shift our data points by an amount -b, 
taking b from the best fits Eq. J29I and Eq. J30t . and use 
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L cc exp(-x 2 (M)/2) (with normalisation appropriately fixed 
to unity) as the likelihood function in Eq. i34\ . where 



X 2 (M) 



s 



Ati - b s 



shift 



a(M)Ki 



CTi 



(35) 



is calculated for all the possible values of a(M), defined by the 
coefficients of K\ (K q ) in Eq. Jl 3i and Eq. dl4> . respectively. 
The sum in Eq. (I35> is taken over the all the data points Af;, 
which are symbolized by T>, with <x,- characterizing the uncer- 
tainties in the measured time lags. The calculated ^ 2 (M L ) for 
the different combinations of the data sets we use in the case 
of linear corrections to the refractive index Eq. fl!4i is shown 
in Fig. [3] The minima of x 2 correspond the 'signal-like' re- 
gions, where the data are better described by a scenario with a 
refractive index, induced by quantum gravity. The most robust 
estimation on the lower limit of quantum gravity with a linearly 
energy-dependent correction is obtained from the combination 
of all the data sets, and is indicated by a solid line in Fig.[S] 



M L > 6.9 • 10 15 GeV 



(36) 



Similar considerations lead to the following lower limit on 
quadratic quantum-gravity corrections Eq. (0 to the photon 
dispersion relation: 



M Q > 2.9 • 10 6 GeV. 



(37) 



To the accuracy stated, we find identical numerical results, 
whether we use a logarithmic measure for M cut off at Al = 
10 18 ' 19 ' 20 GeV, as shown in Fig. [3] or a 1 /M measure integrated 
to infinity, 



I IM 

_0 

oo 



= 0.95, 



(38) 



where L'{l/M) is the likelihood function with respective to 
1/M, as shown in Fig.|9] 

The key result (1361 is significantly stronger than that in 
(Ellis et al. 2000b), thanks to the improved analysis technique 
using wavelets and the use of a more complete dataset. 

8. Discussion 

We have investigated in this paper possible non-trivial prop- 
erties of the vacuum induced by quantum gravity, by probing 
modifications of the dispersion relation for photons. These fea- 
tures can appear in several approaches to quantum gravity, in- 
cluding Liouville string theory, models with large extra dime- 
sions (Campbell-Smith et al. 1999) and loop gravity. Similar 
modifications of dispersion relations can take place for the 
other matter particles, leading to other non-trivial effects such 
as changes in the thresholds for some reaction attenuating ultra- 
high energy cosmic rays (UHECR), or vacuum Cerenkov-like 
radiation (see (Sarkar 2002) and references therein), which 
could have a large influence on the interpretation of the puz- 
zling astrophysical data on UHECR. 



We have attempted to extract the most complete informa- 
tion about the possible vacuum refractive index induced by 
quantum gravity, by using wavelets to look for any correla- 
tion with redshift of the time lags between the arrival times 
of sharp transients in GRB light curves observed in high- and 
low-energy spectral bands. This analysis combined continu- 
ous wavelet transforms to remove noise and discrete wavelet 
transforms to identify sharp transients in different spectral 
bands. Eight GRBs with known cosmological redshifts and 
light curves available publicly have been used in our analysis. 

It is instructive to compare the time lags we find with those 
found using cross-correlation analysis (Band 1997; Norris et 
al. 2000; Norris 2002). Our measurements of the spectral time 
lags in BATSE 64 ms light curves for five GRBs, which are 
common for the sample we used and that under consideration 
in Norris et al. 2000, are consistent in absolute value with the 
trend found in Norris et al. (2000) for bursts with higher lu- 
minosities (which are closer, on the average) to have shorter 
time lags. Discrepancies are found in two cases out of five 
common GRBs. Namely, we found soft-to-hard evolution for 
GRB971214 and GRB990123 (positive time lag), which is op- 
posite to the hard-to-soft evolution (negative time lag) found in 
Norris et al. 2000. These two GRBs havea complicated struc- 
ture of emission: GRB971214 has a wide "clump" of emis- 
sion which consists of spikes that are barely overlapped, while 
GRB990123 consists; in two intensive wide pulses with a quite 
complicated cluster afterwards. Thus these two GRBs can be 
attributed to a lack of morphological classification power of 
the cross-correlation technique due to the problem of interpre- 
tation of the CCF pike width (Band 1997). As an explicit ex- 
ample, we analysed GRB941119, which has been assigned, in 
Band 1987, as the GRB without clear spectral evolution with 
the respective to cross-correlation analyses. The light curve 
of GRB941119 consists in cluster with several closely spaced 
spikes protruding from a smooth envelope. We found five pairs 
of identical genuine variation points in first and third spectral 
bands, demonstrating all together the soft-to-hard spectral evo- 
lution with weghted average time lag Af = 0.051 ± 0.384 s. 
One can see from Table [2 that in several cases different pairs 
of identical genuine variation points detected in one and the 
same GRB's light curve demonstrate different kinds of spec- 
tral evolution, namely either hard-to-soft or soft-to-hard. This 
fact could be connected with a possible intrinsic spectral evo- 
lution during the burst progress. The example of GRB941119 
is one of such cases, demonstrating unclear spectral evolution 
with the respect to the cross-correlation analysis. The wavelet 
technique we apply classifies more explicitly a variable spectral 
evolution and consequently gives more accurate results for the 
weighted average spectral evolution, as in case of GRB941 1 19. 
The time lags measured between the BATSE 25-55 keV and 
OSSE 3-6 MeV light curves exceed substantially the time lags 
between the first and third BATSE energy bands. Without the 
correction by the ratio Eq. (13 1 i . the BATSE-OSSE time lags 
we find are ^ -2 s, in good agreement with Piro et al.1998, 
where the time lag between bands of GRB970228 with a simi- 
lar energy difference was estimated using BeppoS AX data. The 
wavelet technique is very effective also to deal with transient 
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signals such as TTE data with low signal-to-noise ratios, where 
cross-correlation analyses can fail. 

We believe that there is no regular cosmological evolution 
of the sources we use. This fact is already widely accepted 
in the literature on other high-redshift sources, namely super- 
novae. So there is no effect that can cause any correlation of 
weighted average time lags with redshifts of the sources. Thus 
any correlation of spectral time lags we may find can be at- 
tributed to the effects of propagation. 

Finally it should be stressed that the method we applied 
is very powerful in the sense of analysing signals with strong 
non-linear dynamics behind, as GRBs could be. Of course it 
does not pretend to explain the underlying dynamics and phys- 
ical origin of GRBs, but it gives a hint of the stage at which the 
most variable processes take place and characterize quantita- 
tively the degree of instability accompanying those processes. 
In our case the radiation from genuine variation points is con- 
sidered as the messenger of a fast non-linear dynamics at the 
source. The Lipschitz exponents, which characterize quantita- 
tively the 'degree' of instability of the dynamics, give the in- 
formation of whether photons of different energies have been 
produced at one and the same event at the source. So this gives 
an ideal 'time offset' between energy bands to measure any dif- 
ferences in the speed of light. 

It is widely accepted (see Band et al. 1997; Norris et al. 
2000, and references therein) that the spectral evolution in 
GRBs leads to peaks migrating later in time. These time lags 
are not directly connected with the distance to the source, but 
are correlated with intrinsic properties of GRBs, such as lumi- 
nosity (Norris et al. 2000) or variability (Shaefer et al. 2001). 
The quantum-gravity energy-dependent time delay plays the 
role of a foreground effect of opposite sign to the usual spec- 
tral evolution of GRBs, which increases with distance. Hence 
model-independent information about quantum gravity can be 
extracted only from a statistical analysis of sources with a 
known distance distribution. 

We have not found a significant correlation of the mea- 
sured time lags with the cosmological redshift that would in- 
dicate any deviation of the vacuum refractive index from unity. 
This fact allows us to establish 95 % C.L. lower limits on 
the quantum-gravity scale at the levels 6.9 • 10 15 GeV and 
2.8 • 10 6 GeV for linear and quadratic distortions of the disper- 
sion relation, respectively. However, despite the lack of any sig- 
nificant evidence for a quantum-gravity signal, there is a region 
of the linear scale parameter M where the data are better de- 
scribed by a scenario with a refractive index induced by quan- 
tum gravity. This fact indicates that any increase in the statis- 
tics, especially with higher resolution, would be of the utmost 
interest for exploring the possibility of a quantum-gravity cor- 
relation of spectral time lags with redshift. For this reason, we 
urge the light curves for all GRBs with measured redshifts to 
be made generally available, as is already the case for BATSE 
data. 

It has been observed that, if the dispersion laws for el- 
ementary particles differ from the standard ones, the expan- 
sion of the Universe may result in the gravitational creation 
of pairs of particles and antiparticles with very high energies 
(Starobinsky & Tkachev 2002). The expansion of the Universe 
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(both at present and in the early Universe) gradually redshifts 
Fourier modes of a quantum field, and may transport them from 
the trans-Planckian region of very high momenta to the sub- 
Planckian region where the standard particle interpretation is 
valid. Then, if the WKB condition is violated somewhere in the 
trans-Planckian region, the field modes enter the sub-Planckian 
region in a non-vacuum state containing equal numbers of par- 
ticles and antiparticles. The most restrictive upper limit fol- 
lows from the number of UHECR created at the present epoch 
(Starobinsky & Tkachev 2002). This limit, together with our 
measurements of the vacuum refractive index, may rule out the 
possibility of detecting imprints of trans-Plankian physics on 
the CMB anisotropy, as proposed in (Brandenberger & Martin 
2002). 

It is a widespread belief that the combination of quan- 
tum theory with gravity is an explosive area, to the extent 
that it blows up the concrete fabric of space-time, splitting it 
into space-time foam. It is true that, as yet, there is no uni- 
versally accepted mathematical model of this metamorphosis, 
while several attempts at it are characterized by a varying level 
of success. However, despite the firm and widely-held opin- 
ion that quantum gravity defies experimental tests, it is re- 
markable that several such tests have been proposed in recent 
years. Some of them have put rather severe constraints on some 
quantum-gravity models, whilst some models based on non- 
critical strings seem to remain unscathed so far (Ellis et al. 
2002; Ellis et al. 2002a). In particular, as we have shown in this 
paper, using the most sophisticated technique available, that of 
wavelets, to pick up authentic time-lag effects in gamma-ray 
bursts, we have been able to put a rather firm lower bound on 
the dispersion of light in the vacuum: M > 6.9 x 10 15 GeV. This 
approaches the range of scales where we might expect such an 
effect to appear. Data able to test further this possibility already 
exist, and more data will soon be available. 

Une affaire a suivre .... 

Acknowledgements. This work is supported in part by the European 
Union through contract HPRN-CT-2000-00152. The work of D.V.N, 
is supported by D.O.E. grant DE-FG03-95-ER-40917. We are grate- 
ful to A. Biland for his kind help for reading out the TTE data. We 
thank Hans Hofer for his continuing support and interest. We are also 
greateful to I.M. Dremin for useful conversations on wavelets, and 
A.S.S. is indebted to E.K.G. Sarkisyan for many useful conversations 
on statistics. 

Appendix A: Wavelet-Based Thresholding 
Estimator 

A thresholding estimator of a n = 2 '-component discrete signal 
X[n] decomposed at the resolution level L (see Section 4) in a 
wavelet basis can be written (Mallat 1998) as: 

.1-1 2'-' 2 L -' 

P = Z Z PT0^pn + J] PTrfrnXPlm, (A. 1) 

j=L m=0 m=0 

where the function p T provides a soft threshold: 

( x - T if x > T 
p T (x) = I x + T if x < -T (A.2) 
(0 if \x\ < T 
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If the signal X[n] is contaminated by additional noise W[n], the 
threshold T is generally chosen so that there is a high proba- 
bility that it is just above the maximal level of this noise. Since 
W[n] is a vector of N independent Gaussian random variables 
of variance cr 2 , one can prove (Mallat 1998) that the maximum 
amplitude of the noise has a very high probability of being just 
below T — cr -J2 logn. 

A signal X[n] of size n has n/2 wavelet coefficients 
{dj-i, m }o< m <n/2 at the finest scale. The coefficients of / itself, 
\dj m \, in the sum Eq. \2\\ are small if / is smooth over the sup- 
port of <£r/-i, m , implying dj_ lm ~ dj_ lm . In contrast, \d f } _ l m \ 
is large if / has a sharp transition in the support of \j/j-i,m- A 
piecewise-regular signal has few sharp transitions, and hence 
produces a number of large coefficients that is small compared 
to n/2. At the finest scale, the signal / thus influences the value 
of a small portion of large-amplitude coefficients dj-i, m , that 
are considered to be 'outliers' . All the others are approximately 
equal to dj_ x , which are independent Gaussian random vari- 
ables of variance cr 2 . 

A robust estimator of cr 2 is calculated from the median of 
(dj X )o<m< n /2- The median of P coefficients Med(a p )o< p< p is 
the value of the middle coefficient a„ of rank P/2. As opposed 
to an average, it does not depend on the specific value of coef- 
ficients a p > a m . If M is the median of the absolute value of P 
independent Gaussian random variables of zero mean and vari- 
ance cr 2 , then one can show (Mallat 1998) that the expected 
value of M is 0.6745ctq. Thus the variance cr 2 of the noise W 
is found from the median Mx of (\dj ^ J)o<m<«/2 by neglecting 
the influence of /: 

cr=J^. (A.3) 
0.6745 V 

One may say that / is responsible for few large amplitude out- 
liers, and that these have little impact on Mx- In practice, it is 
more convenient to transform the signal, by scaling it in such 
a way that the wavelet coefficients at the finest level of decom- 
position have median absolute deviation equal to unity, as seen 
in Fig |2 . 

Appendix B: Lipschitz Regularity 

Taylor expansion relates the differentiability of a signal to local 
polynomial approximations. Suppose that / is m times differ- 
entiable in [v — h; v + h]. Let p v be the Taylor polynomial in 
the neighborhood of v. The m order of differentiability of / 
in the neighborhood of v yields an upper bound on the error 
£y(0 = f(t) - Pv(t) when t tends to v. The Lipschitz regularity 
extends this upper bound to non-integer exponents. Namely, a 
function / is pointwise Lipschitz: a > at v, if there exists a 
polynomial p v of degree at most a such that 



continuosly differentiable in this neighborhood. If v < a < 1, 
then p v (t) = /(v), and the Lipschitz condition Eq. JB. 1 1 be- 
comes 



\f(t)-p v (t)\<'K\t-y\ a 



(B.l) 



where 7C is a constant. The polynomial p v (t) is uniquely de- 
fined at each v. If / is m < a times continuously differentiable 
in a negborhood of v, then p v is the Taylor expansion of / at 
v. Pointwise Lipschitz exponents can vary arbitrarily from ab- 
scissa to abscissa. Iff is uniform Lipschitz: a > m in the neigh- 
borhood of v, then one can verify that / is necessarily m times 



1/(0- /(v)|<7C|*~vr, 



(B.2) 



/ is not differentiable at v and a characterizes the singularity 
type. For more references to the mathematical literature, see 
(Dremin et al. 2001; Mallat 1998). 
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GRB (BATSE Trigger) 


z 


Time lag, error (s) 


970508 (6225) 


0.835 


BATSE (64 ms) 

Af BATSE = _ .064 ±0.860 


971214 (6533) 


3.418 


BATSE (64 ms) 
Ar BATSE = o.i3 +0.860; Af BATSE = -0.064 ±0.860 
cobmbined BATSE (64 ms) 

a? batse _ 033 +0.608 

comb 


980329 (6665) 


3.9 


BATSE (64ms) 

Af BATSE _ o.o64 ±0.860; A^ ATSE = ±0.860 
combined BATSE (64 ms) 

Af c<^b E = °- 032 ±0-608 
BATSE-TTE (2.7 ms) 
A/™ = -0.34 ±0.019 
OSSE (64 ms) rescaled 
Ar ossE _ _o.048 ±0.019; Af 0SSE = ±0.038 
combined BATSE+OSSE+TTE 

A Cmb E+rm0SSE+TTE = -°- 036 ±0 - 013 


980425 (6707) 


0.0085 


BATSE (64 ms) 

Ar BATSE = _j 792 + L 705; A^ ATSE = -1.28 ±0.860 
combined BATSE (64 ms) 
A db E = -1-384 ±0.768 


980703 (6891) 


0.966 


BATSE (64 ms) 

a? batse _ _ .832 ±0.860 
OSSE (64 ms) rescaled 
Ar ossE = _ .040 ±0.019 

combined BATSE+OSSE 

Af BATSE+OSSE = _ Q m ±Q Ql 9 
comb 


990123 (7343) 


1.600 


BATSE (64 ms) 

ArfATSE _ o.230 ±0.860; Af BATSE = -0.064 ±0.860; Af BATSE = -0.128 ±0.860 
combined BATSE (64 ms) 
Af BAT b SE = 0.013 ±0.496 
OSSE (64 ms) rescaled 
A; osse = _ .049 ±0.019; Af 0SSE = -0.046 ±0.019; Af° SSE = -0.045 ±0.019 
combined BATSE+OSSE 
A S E+ ° SSE = -°- 047 ±0 - 011 


990308 (7457) 


1.2 


BATSE (64 ms) 

Af BATSE _ o +0.860; Af BATSE = -0.064 ±0.860; A^ ATSE = -0.256 ±0.860 
Aj batse = _j 024 ±0.860; At% ATSE = ±0.860; A^ ATSE = 0.064 ±0.860 
Af BATSE = o +0.860 
combined BATSE (64 ms) 

A CTb SE = -°- 183 ±°- 325 


990510 (7560) 


1.619 


BATSE (64 ms) 

AfB ATSE = 0.384 +0.860; AfB ATSE = 0.448 ±0.860; Af BATSE = ±0.860 
A( batse _ _ .256 ±0.860; Af BATSE = ±0.860; Af BATSE = -0.528 ±0.860 
Af BATSE _ _ .128 ±0.860; Af BATSE = -0.256 ±0.860; Ar BATSE = ±0.860 
AfBATSE _ _ .448 +0.860 
combined BATSE (64 ms) 

A CST = -°- 078 ±0 - 272 

OSSE (64 ms) rescaled 

Af ossE _ _ .045 ±0.019; Ar° SSE = -0.032 ±0.019; Ar° SSE = -0.041 ±0.019 
combined BATSE (64 ms) 

AfBATSE+OSSE = _q m g ±Q Qn 
comb 


991216 (7906) 


1.02 


BATSE (64 ms) 

Arf ATSE = -0.064 ±0.860; A<f ATSB = -0.064 ±0.860; A<f ArsB = -0.064 ±0.860 
Af BATSE = -0.064 ±0.860; At% ATSE = 0.064 ±0.860; Ar BATSE = ±0.860 
Ar BATSE = o +0.860; Af BATSE = ±0.860 
combined BATSE (64 ms) 
AfBATSE _ _o.040 ±0.304 

comb 



Table 1. Data on the light curves for GRBs with known redshifts used in this analysis. The third column gives the time lags 
between the arrivals of every identical genuine variation points in the third (high-energy) spectral band and the first (low- 
energy) one. For every light curve, we give the weighted means of the time lags in 64 ms BATSE domain, combining the genuine 
variation points for that GRB. Both statistical and systematic errors are included. The average spread of individual time lags is 
below the mesurement uncertainty. Also indicated are the results obtained by combining OSSE 64 ms light curves in the 3-6 MeV 
energy range with BATSE light curves in the 1 15-320 keV energy band. In the case of GRB980329, the results of BATSE 64 ms 
and TTE 2.7 ms resolution measurements are combined with the OSSE-BATSE6A ms comparisons. 



